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Abstract 


This investigation is concerned with the dynamic synthesis of 
gyroscopic structures consisting of point-connected substructures. 

The object is to develop a mathematical model capable of an adequate 
simulation of the modal characteristics of a helicopter using 
a minimum number of degrees of freedom. The basic approach is to 
regard the helicopter structure as an assemblage of flexible sub- 
structures. This approach has the advantage that it permits the 
representation of each substructure by a given number of "modes". 

The variational equations fox the perturbed motion about certain 
equilibrium solutions are derived. Assuming that some of the sub- 
structures are continuous, the variational equations are of the hybrid 
type so that it is necessary to discretize them, a task done quite 
well in the context of the component-mode synthesis. The discretized 
variational equations can be conveniently exhibited in matrix form, 
and a great deal of information about the system modal characteristics 
can be extracted from the coefficient matrices. The derivation of the 
variational equations requires a monumental amount of algebraic opera- 
tions. To automate this task a symbolic manipulation program on a 
digital computer is developed. 
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1. Introduction 


This investigation is concerned with the dynamic synthesis of 
helicopters. The object is to develop a mathematical model capable of 
an adequate simulation of the modal characteristics of a helicopter 
using a minimum number of degrees of freedom. 

The method of approach represents an extension of the component- 
mode synthesis (Refs. 1 and 2) to gyroscopic systems. The basic idea is 
to regard the helicopter structure as an assemblage of flexible substruc- 
tures. In the case of a helicopter, the various substructures can be 
identified as the airframe, the transmission shaft, the main rotor, and 
the tail rotor. The rotor can be slmul&ted at a single substructure or 
as a collection of flexible blades. To ensure that the various substruc- 
tures are acting as a single system, an orderly klnematical procedure 
for point-connected substructures is developed (Ref. 3). 

The substructure synthesis approach has the advantage that it 
permits the representation of each substructure by a given number of 
"modes". This number can be varied according to the stiffness properties 
of the substructure. Clearly, a stiff substructure can be represented 
by a smaller number of degrees of freedom than a flexible one. This is 
not the only advantage, however. Indeed, the approach has the added 
versatility that it permits a given substructure to be represented by a 
discrete mathematical model or a continuous one. For example, it is 
reasonable to represent the fuselage by a discrete model and a rotor 
blade by a continuous model. Assuming that some of the sub structure *7 
are continuous, the system equations of motion are hybrid in the sense 
that they comprise both ordinary and partial differential equations. 

Note that even in the absence of discrete substructures the system in- 
cludes ordinary differential equations describing rotational motion. 

The equations of motion can be conveniently obtained by the Lagrangian 
approach. The general equations are nonlinear and solutions of hybrid 
systems of nonlinear differential equations are beyond the state of the 
art. Fortunately, general solutions are not really necessary, nor are 
they likely to be very informative. Indeed, quite often one is content 
with solutions in the neighborhood of given special solutions, where the 
latter are often called equilibrium solutions. The motion in the neigh- 
borhood of the equilibrium solutions is referred to as the perturbed 
motion. When the perturbed motion is small, the linearized equations 
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obtained by ignoring the higher-order terms are known as the variational 
equations (Ref. 4). Two special solutions of considerable Interest In 
Che area of helicopter dynamics are associated with hover and forward 
flight. 

The variational equations are still of the hybrid type. For a 
computational treatment of the equations, it is necessary to discretize 
them. This can be done quite well in the context of the component-mode 
synthesis. The procedure amounts to representing a continuous substruc- 
ture by a finite series of space-dependent modes multiplied i>y time- 
dependent generalized coordinates. The meaning of these modes repre- 
sents a study in itself and there appears to be some measure of freedom 
in their selection. The discretized variational equations of motion can 
be conveniently exhibited in matrix form, and a great deal of informa- 
tion about the system modal characteristics can be extracted from the 
coefficient matrices. 

In the case in which one of the substructures rotates uniformly in 
the equilibrium state, such as in the case of the helicopter rotor, the 
coefficient matrices of the variational equations are generally periodic 
functions of time. However, when the rotor possesses both Inertial and 
structural symmetry, which is possible if the rotor has at least three 
blades, a coordinate transformation can reduce the periodic matrices of 
the variational equations to constant ones. The procedure is sometimes 
referred to as the method of multiblade coordinates (Ref. 5). The 
resulting constant-coefficient equations represent a typical linear 
gyroscopic system. 

The derivation of the variational equations of motion and their re- 
duction to a constant-coefficient system requires a monumental amount of 
algebraic operations and differentiations. In addition, there is the 
question of producing an ordering scheme to derive equations cf consis- 
tent order (Ref. 6). This points to the necessity of automating the 
various derivations, a task accomplished by means of a symbolic manipula- 
tion program on a digital com’-i ter . 

In this report, the above aspects are first treated in a somewhat 
general way and then applied to a simplified mathematical model of a 
helicopter . 
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2. Kinematics! Considerate ns 

Let us consider a general structure consisting of a central sub- 
structure C and a given number of appended substructures (see Fig. 1). 

The appended substructures can be of many types. However, fur the pur- 
pose of this paper we shall confine our discussion to a few specific 
types, namely, flexible substructures rotating relative to the central 
substructure (type A) and chains of substructures composed of a flexible 
substructure rotating relative to the central substructure to which 
another flexible substructure B is connected (type AB) . The various types 
are illustrated in the illustrative example of Sec. 9. Clearly, rigid 
substructures are simply special cases of flexible substructures obtained 
by ignoring the flexibility. In addition, although we recognize that 
there may be more than one appendage of each type, we shall confine our 
discussion to a representative one of each type, with summation implied 
over the entire number of substructures of the same type. The kinemati- 
cal relationships pertinent to such general structures are presented in 
detail in Ref. 3 and are only reviewed briefly here. 

Let us consider the inertial system of axes XYZ with the origin at 
0 and identify a system of axes x^y^z^ with the origin at an arbitrary 
point C of the central substructure. Then, denoting the velocity vector 
of point C relative to point 0 in terms of components along XYZ by 


w QC> the angular velocity of the frame x c y c z c relative to XYZ in terms 
of components along x c y c z c by and the matrix of direction cosines 
between x c y c z c and XYZ by T CQ , the absolute velocity of any mass point 


in the central substructure in terms of components along x n y n z r is 

"c = T co~oc + n c ( ^c + + “c = T co?oc ' (t c + u c>?c + “c (1) 


In Eq. (1) r^, is the position vector and u^, is the elastic displacement 
vector of the mass point measured relative to x y z ., Q is a skew sym- 
metric matrix associated with ft (see Ref. 4), ? + u is a skew sym- 

^ v L 

metric matrix associated with r c -f u c> and u c is the elastic velocity of 
the point relative to axes x y z . 

V V V 

Next, we consider the appended substructure A attached to the cen- 


tral substructure at point A and rotating relative to it. First, we 


introduce a nonrotating set of axes x CA y CA z CA in substructure C with 
origin at the point A and moving relative to axes x y z as substructure 
C deforms. Then, introducing another set of axes x.y z attached to 

AAA 
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substructure A with origin at A and moving with substructure A relative 
to substructure C, we can write the absolute velocity of any mass point 
In substructure A in terms of its elastic velocity relative to the frame 


x.y.z., the motion of the frame x.y.z. relative to axes x_.y_.z_., and 
AAA AAA CA CA CA 


the absolute motion of axes x„. y„ A z„. . To this end, we recognize that 

uA LA LA 

the absolute velocity of the attachment point A is obtained simply by 
evaluating Eq. (1) at the specific point A. We denote this velocity by 
w^. The absolute angular velocity of axes x cA y CA Z CA 

( 2 ) 


Oca ~ l ea l gaOc + ~ca 


where is the matrix of direction cosines between axes x c ^y C A z CA and 


X C^C Z C ^ Ue tQ £ eometr y alone, is the matrix of direction cosines 


between these axes due to the deformation of substructure C, and 3 


CA 


is the angular velocity of X £^y C A z cA relative to x^y^,z c due to the 


deformation of substructure C. Explicitly, 3^ is the time derivative 
of the linear rotation vector 

(3) 


?CA L GA ( ' V C~CA ) 


of axes x y z relative to axes x y z in which V is the skew sym- 
LA LA LA L L L L* 


metric operator matrix corresponding to the curl operator associated 


with the rectangular coordinate system x y z ; the subscript A of u 

L L L ~CA 


specifies that the quantity in parenthesis is to be evaluated at point 


A. In addition, the matrix L„. has the form 

LA 

L EA 1 + g CA 


(4) 


where 1 is the unit matrix and 3 CA is the skew symmetric matrix associated 
with the vector $ . Denoting by ui. the angular velocity of body axes 

~LA ~ A 

x^y^z^ relative to axes ^^CA^A an d letting be the matrix of direc- 
tion cosines between these axes, the absolute angular velocity of axes 
x^y^z^ in terms of components along these axes can be written in the 


form 


?A L aC‘CA + ~A T AC?C + L A~CA + “A 


(5) 


where 


AC 


l a l ea l ga 


( 6 ) 


is the matrix of direction cosines between axes x^y^z^ and x^y^z^. 


Finally, by analogy with Eq. (1) the absolute velocity of an arbitrary 
point in A is 


~A T AC~CA + ^A^A + ~A^ + “A T AC~CA “ (r A + U A J ?A + ~A (7) 






1 1 . i L 1 - I 


Turning our attention to the specific chain of substructures, we 
recognize that the first substructure in the chain is exactly analogous 
to a type A substructure. To describe the motion of the second substruc- 
ture, we first consider the reference frame x.„y.„z.„ in substructure A 
with origin at the point of connection, point B, between the second sub- 
structure B and the first substructure A and which moves relative to the 
frame x A y A z A as substructure A deforms (see Fig, 1). Next, introducing 
the axes x^y^z^ with the origin at point B and attached to substructure 
B, we can write the absolute velocity of any mass point in substructure 
B in terms of its velocity relative to the frame x„y z , the motion of 
the frame x D y D z D relative to x.^y.^z. and the absolute motion of the 
frame By analogy with Eq. (5), we can write the absolute 

angular velocity of axes x n y„z in terms of components along these axes 

D D D 

in the form 


?B T BA?A + L B?AB + _B 


( 8 ) 


where ft A is given by Eq. (5), is the angular velocity of axes x gy 0 z B 

relative to x A 0 y AB z AB » T BA = L b L EB L GB is dhe matrix of direction cosines 

between axes x gy B z B and x A y A z A in which Lgg is due to the geometry and 

L is due to the deformation of body A, and by analogy with Eq. (3). 

Ed 


!ab ~ l gb ( 7 b“ab^ 


(9) 


is the linear rotation vector of axes x A 0 y AB z A g relative to x A y A z A - 
Moreover, by analogy with Eq. (7), the absolute velocity of an arbitrary 
point in B is 


?B * T BA”AB + V^B + V + "B ' T BA"AB * <r B + “b^B + ?B 


( 10 ) 


where w AD is obtained by evaluating Eq. (7) at the point B and u is 
^ Ad - d 

the elastic velocity relative to axes x^y^z^. 


3 . System Lagrange's Equations in General Hybrid Form 

The object is to derive the system Lagrange's equations of motion 
in general form. To this end, it is necessary to determine the func- 
tional dependence of the Lagrangian. If summation over all substruc- 
tures is implied, then we can write the system kinetic energy in the 
general form 


T 


1 

2 



•T* ^ 

”c?c dm c 




•T* J 
w.w. dm. 
_A_A A 



•I' a 
?B”B dm B 


(ID 
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where v c , w^, and w 0 are given by Eqs. (1), (7), and (10), respectively. 

Assuming that in equilibrium substructure A and substructure B rotate 

with the uniform angular velocities id. about z. and about z_, respec- 

A A O D 

tively, while any other motion is zero, we can write 


5fc • e c- e c 




“a5a + 0 A?A 


( 12 ) 


5« “ “b^b + 0 B?b 


where is the vector of direction cosines between z^ in the equili- 
brium configuration and axes x A y A z A a nd 2.^ is the vector of direction 
cosines between z_ in the equilibrium configuration and axes x_y 0 z_. 

a ODD 

Moreover, 0 C , 0^, and 0 & are 3*3 matrices depending on the oscillation 
of axes x c y c z c relative, to XYZ, x A y A z A relative to x^y^z^, and 
relative to x A gy A £ Z A g » respectively, and 0^, 0^, and § B are vectors of 
generalized velocities. Hence, the functional dependence of the kinetic 
energy is 

T = T( ”0C’ 9 Ci’ ic’ "C’ “c* ~CA’ ~CA’ 6 Ai’ ®A» ~A’ V HaB’ “aB» 


3 Bi’ ?b* i?B* ^ 


(13) 


where u„. denotes u„ evaluated at point A of substructure C and u._ 
denotes u^ evaluated at point B of substructure A. 

Assuming that the potential energy is due entirely to elastic 
effects (including possible elastic restraints between substructures), 
we shall write it in the general form / 

V = U(9 , 


Al* A2 


t-u Wiiv- gctiCLax turui - 

..2’ 0 B1’ 6 B2 ) + 2 X c ~C A DC-C dD C + ~C A SC~ 

+ 2 Jd ~A A DA~A dD A + 2 ~A A SA~A dS A + 1 jv “b^DB-B 

A A B 

+ 2 X ~B A SB~B dS B 


where U represents the potential energy due to elastic restraints of the 
rotation of substructure A relative to C and substructure B relative to 
A and A^, A gc , A^, A g ^, A QB and A gB are "two-sided" differential 
operators, containing partial derivatives with respect to the spatial 
variables. The operators are symmetric and positive definite and can be 
identified as energy operators over the respective domains. The last 
six terms of Eq. (14) represent symbolic operations involving integrations 


over the elastic domains D, , and D_ and their boundaries S , S., and 

LAB LA 

Sg, respectively. 

Next, introduce the Lagrangian l = T-V and write it in the form 


L = 


L c v d c + X A v dD A 

*Js K i SB dS B - U 

D 



L DB dD B 




L sc ds c 



L SA dS 


A 


(15) 


where 

’ L sc’ 

l da* 

l sa* l 

Hamilton* i 

s principle 

(Ref. 

obtain Lagrange 

* s equations 

— / 

3L \ 

= 0 


dt \ 

3w„„ ) 
-0C / 

* 

3L 

_ <L | 

f 3L 

)-°. 

Cx 

dt < 

Ki. 




3L ( 



36 Bi ' 

3£ dc 


'*hc 

) + L c u . 

9 “c 

3t \ 

- 

s£ da 

3 "a 

- -I 

at ! 

^Da' 
^?A . 

) + 

^DB 

_i_/ 

' 9£ db' 

) + _ l b“i 

3 “b 

at \ 

k 3 _b 


"SB 



3L 

d 

( 9L V 

J 

"Ai 

dt 


d 


1= o 

i=l,2,3 

dt 

\ 39 bJ 



over D„ 


(16a) 


(16b) 


(16c) 


where L , L , and L are linear differential operator vectors reflect- 
^ A ^ 15 

ing stiffness properties. Equations (16c) are subject to appropriate 
boundary conditions, which are omitted for brevity. Note that the 
relation between the operators A and L is given by 

(17) 


J ^ u T A Q udD + J s u^AgudS = u 1 
where the satisfaction of the boundary conditions is implied. 
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4. Variational Equations 

The hybrid set of equations, Eqs. (16), is generally nonlinear, 
particularly in the rotational coordinates. Moreover, due to the uni- 
form angular velocity tu A , the angle to A t appears in the matrix of direc- 
tion cosines giving the orientation of axes x ^y^ z ^» a nd a similar state- 
ment can be made for axes x^ygZg. Therefore, Eqs. (16) possess time- 
dependent coefficients, or more specifically periodic coefficients. 
General solutions of sets of hybrid nonlinear differential equations 
with periodic coefficients do not exist. However, our interest is in 
solutions in the neighborhood of known solutions. These known solutions 
are either constant solutions or periodic solutions and define an equili- 
brium state or equilibrium motion, respectively. To study the dynamic 
characteristics in the neighborhood of known solutions, the perturbation 
equations are derived. The linearized perturbation equations are known 
as variational equations. To derive the variational equations, one can 
expand the Lagrangian in a Taylor series. The linear terms define the 
equilibrium and can be ignored since they cancel out in the equations of 
motion. The quadratic terms lead to the linear terms in the equations 
of motion and are the only terms retained in the Lagrangian. Note that 
this procedure involves simply a transformation of coordinates so that 
the equilibrium of the new equations is trivial. However, the Lagrangian 
obtained in this way will still possess periodic coefficients. This 
creates a host of problems, as far as the analysis of the variational 
equations is concerned. Indeed, systems of equations with periodic 
coefficients are considerably more difficult to analyze and the results 
considerably harder to interpret compared with systems with constant 
coefficients. Under certain circumstances, it is possible to reduce a 
system with periodic coefficients in time to one with constant coeffi- 
cients. Such systems are known as "reducible" and the reduction process 
involves a linear transformation, where the equations of the transforma- 
tion also have periodic coefficients with the same period as the time- 
dependent terms to be eliminated. However, there is no general analyti- 
cal procedure for determining the linear transformation leading to the 
constant-coefficients system. We shall not be concerned with a general 
reducible system but with a special kind, namely, one for which a real 
transformation exists such that the transformed Lagrangian remains real 
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and does not depend explicitly on time. Transformations with these 
properties can easily be found if the rotating substru turea A and B 
possess certain symmetry. For example, if substructure B consists of 
three or more Identical beams with one end fixed to a hub and spaced 
equally around the hub in a plane perpendicular to the axis of rotation, 
and substructure A consists of a uniform shaft rigidly attached to the 
hub, then the transformation can be found by the method of multiblade 
coordinates (Ref. 5). In the following, we shall assume that transfor- 
mations u* = u*(u A ,t), u* = u*(u_,t) exist such that the transformed 
Lagrangian does not depend explicitly on time. Taking into account the 
above discussion, the linear constant-coefficients equations of motion 
can be written in a form resembling Eqs. (16). In fact, the form remains 
the same if the generalized coordinates entering into Eqs. (16) can be 
regarded as measuring displacements from equilibrium u 

Next, we must recognize that a system of the type under considera- 
tion possesses Ignorable coordinates. Our object being to study the 
modal characteristics of the system, it is desirable to eliminate these 
ignorable coordinates. To this end, the three components of translation 
are r.-cogiized immediately as ignorable coordinates. It is not Immediately 
apparent, however, that certain rotational coordinates may also be 
ignorable. Indeed, if we assume that in equilibrium the z c axis is 
aligned with the z^ axis, then the angles 0^ and 6 ^ , defined as the 

rotation of substructure C about the z axis and the deviation of the 

L» 

rotation of substructure A about the z^ axis from steady spin, respec- 
tively, are recognized immediately as being ignorable coordinates. 

Since the generalized momenta associated with ignorable coordinates are 
conserved, we can regard the conservation of these generalized momenta 
as constraints on the system and use these constraints to eliminate the 


vector Wq£ and the rotational coordinates 0^ a °d 0^ from the formula- 
tion. A similar discussion leads to the elimination of 0^* 1° view of 

these constraints, the final variational equations can be written in a 


form resembling Eqs. (16b) and (16c) if in Eqs. (16b) we delete i=3. 


5. System Discretization 

The system of equations (16b) and (16c) is hybrid and not very 


convenient for studying its dynamic characteristics. The reason can be 
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easily traced to the fact that some of the coordinates depend not only 
on time but also on the spatial position. Hence, we shall find it 
necessary to eliminate the spatial dependence. This we shall do by a 
process referred to as series discretization. Before proceeding with 
the discretization, we shall present certain mathematical preliminaries. 

Let us consider a function u that is continuous in the closed 
domain D=D+S, where D denotes the open domain and S Its boundary. The 
function Is assumed to be quadratically integrable in the Lebesgue 
sense, which Implies that the integral 

|fu|| 2 = (u, u) »/u 2 dD (18) 


exists, where ||u|| is known as the norm of u. Our interest lies in 
approximating the function u by some other function. To this end, we 
consider the sequence of ?' notions u^, ...» u^, ... defined over D. 

Then, if lim | |u^ - u| | » 0, it Js said that conv erges in the mean 


Next, let us assume that the function u can be represented in the 


mean, 


domain D by the series u = a j^j’ * s convergent in the 

and consider the partial sum 

U N - £ (W) 


Then, the set of functions is said to be complete if it is possible 
to find an integer N and a set of coefficients a^, a 2 , .... a^ such that 
for any e>0, | ju - UjJ | < e. The question remains as to the nature of 
the functions <}k 0=1,2,,...). To answer this question, we wish to delve 
a little deeper into the properties of the differential operators A and 
L introduced earlier in conjunction with vector functions. In the first 
place, we should point out that, whereas L is an actual differential 
operator, A is an operator only in a symbolic sense. If L is a self- 
adjoint operator, then for any two functions in the field of definition 


of L 

(i-U,v) 



vludD = 



ulvdD = (u, Lv) 


( 20 ) 


Integrating by parts, and considering the boundary conditions of the 
problem, we can write (see analogy with Eq. (17)) 
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( 21 ) 


fo vLudD m fv vAjjUdD + 


s vA s udS 


where and Ag are "two-sided" symbolic operators symmetric In u and v. 
For convenience, we shall use the notation 

m f» vA D udD + X 


[u, v]. 


r s vA s udS 


(22) 


The operators L and A are said to be positive definite if the inequality 
(u, u] A > 0 is satisfied. 

In our particular case, [u, u]^ defines twice the potential energy 
associated with a given, closed domain D. Then, the quantity Ju| ^ 
defined by 


M 


A = lu ’ U] A 


(23) 


is known as the energy norm of u. When it is clear that the norm 
involves the operator A, the subscript A can be ignored. The sequence 
of functions u^, U£» u^, ... converges in energy to u if 

lim - u | =0. If u^ represents the partial sum (19), then if for 


N-x» 


any e>0 it is possible to find an integer N and a set of coefficients 
a^ such that |u - u^| < c, the set of functions <)>^, <J> 2 » 
is said to be complete in energy . 


a ^» a 2» • • • » 


• • * » * 

To generalize these ideas, let L be a positive definite operator 
defined on some Hilbert space H with the scalar product (u, v) and the 
norm defined by Eq. (18). But the scalar product (22) also defines a 
Hilbert space. The Hilbert space defined by the scalar product (22) is 
called complete if any sequence of its elements satisfying the condi- 
tion lim K- ♦„!- 0 has a limit (e.g. is the limit of the sequence 

k , n-*«° 

<f>. if lim I* “ I = 0); otherwise it is called incomplete. In general, 
it is incomplete and must be made complete by adding certain new elements 


to it where the added elements are called limiting elements and the 
limiting element <j> is defined by the sequence $ ^ . The scalar product of 
the limiting elements $ and iJj is defined by 


[4» , = lim [<*> 


n-**> 




and the norm of the limiting element $ by 

!$! = ii® 1 6 | 


n -K» 


(24) 


(25) 
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The new Hilbert space thus constructed is called the energy space and is 
denoted by H^. Note that Eqs. (24) and (25) extend the definition of 
scalar product and energy norm to include functions which are not in the 
field of definition of the operator L. Moreover, because L is positive 
definite, all elements in the space also belong to the initial Hil- 
bert space H and convergence of a sequence in energy implies that it 
also converges in the mean. Hence, the functions 4^ in the partial sum 
(19) belong to the space H^. The question is as to what distinguishes 
the functions belonging to H A from those in the field of definition of 
L. It is clear that the integration by parts leading to Eq. (21) lowers 
the order of the operator involved by a factor of two. Hence, if L is 
of order 2p, then A^ is of order p, so that the requirements on the 
differentiability of the functions 4j are lowered accordingly. More 
importantly, however, the energy integrals, Eq. (22), Cake automatically 
into account the natural boundary conditions, so that the functions 4j 
need satisfy only the geometric boundary conditions. Such functions are 
sometimes referred to as energy functions . In more familiar terminology, 
we shall refer to functions belonging to H^ as admissible functions and 
those in the field of definition of L as comparison functic is (Ref. 7). 

The above concepts can be extended to a vector function u, by sim- 
ply defining the norm as j|u|| = (u, u) = .f u udD as well as the 
energy norm |u| = [u, u] = J ^ u A^udD +7g u A,,udS where once again 

the integrals represent symbolic operations. 

In view of the above, we shall represent the elastic motions in the 
domains D c , D^, and D g by the series 


N C 


NA 


» c - SjVcj * *C.V “A ■ £ 


= 


?B 


-If 


(26) 


?Bj?Bj $ Bl?B 


where <I> , and 4> are 3*N„, 3*N , and 3*N matrices of space-depen- 

dent admissible functions and n^, n^» and are vectors of time-depen- 

dent generalized coordinates. This permits us to introduce the n-dimen- 

T T T T 

sional configuration vector q = E 0 cl 9 C2 e Al 6 A2 0 Bl 0 B2 n C r| A T1 B^ ’ where 
n = 4+N-.+N which consists of generalised displacements measured from 

L A li 
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( 27 ) 


equilibrium. Then, the system kinetic energy can be written as 
1 tT • T • IT T* T 

T * 2 q Mq + q Fq + y q l^q + R q + S q + constants 

in which M and K^, are symmetric matrices. Recognizing that Eq. (27) 
represents a Taylor series expansion of the kinetic energy T about an 
equilibrium state, we conclude that the last three terms in Eq. (27) 
actually define the equilibrium state, with the inference that R is a 
constant vector and S is the n-dimensional zero vector. In addition, 
the constants in Eq. (27) can be ignored because they do not enter into 
the equations of motion. Hence, the last three terms in Eq. (27) can 
simply be omitted. 

In terms of the system configuration vector q, the system potential 

energy takes the quadratic form 
1 T 

v - j q\s 


(28) 


where 


Ky = block-diag [0 i 


2 

ru 


30 39 

Ai Aj 


, 30 Bi 30 Bj . 


X 


T I 

& A 4> dS 
C SC C C l 


X 


\ W d a 


h'L 


. X 


*C W D C + 


4 a Ia saV s a + 


a d 


B A SB W 


±,i = 1. 2 


(29) 


is a symmetric matrix. Hence, the discretized equations of motion can 
be written in the matrix form 

Mq + (F T - F)q + - K^q = 0 (30) 

where M and K^-K^ are positive definite symmetric matrices and F T -F is 
a skew symmetric matrix. Equation (30) represents a typical gyroscopic 
system. The natural frequencies and natural modes of the complete 
structure and the closed-form solution of Eq. (30) can be obtained by 
methods such as those developed in Refs. 8 and 9. The interest here is 
not so much in the response as in the dynamic characteristics of the 
system, and in particular, the truncation effect on the predicted char- 
acteristics. 

6. Inclusion Principle and Truncat i on Considerations 

Following the procedure of Ref. 8, we introduce the 2n-dimensional 
•T T ,T 

state vector x(t) = [q (t) q (t)j , as well as the 2n * 2n matrices 
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f M i O' 

-T-- ’ 

0 I Ky-K^- 


r f t -p « VVI 

------- 

Ll^-Ky I 0 


where 0 is the null matrix of order n and write the eigenvalue problem 

associated with the system (30) in the real symmetric form 

to 2 Iy = Ky , u) 2 Iz *= Kz r = 1, 2, n (32) 

r £r £r r ~r ~r 

where cj r are the natural frequencies of oscillation of the system and 

y and z are the real eigenvectors, corresponding to the real and 
~ r ~ r X -1 

imaginary parts of x r » Moreover, K = G I G is a real symmetric matrix. 

Assuming that X is positive definite, it follows that K is positive 

definite, so that the eigenvalues are not only real but positive. In 

2 

addition, the eigenvalues (r*l t , 2, ...» n) have multiplicity two and 

because I and K are positive definite all the eigenvectors are indepen- 
dent and can be rendered mutually orthogonal with respect to the matrix 


Next, let us use the Cholesky decomposition and write X in the form 
X 

I^LL , where L is a lower triangular matrix. Introducing the notation 
T T 

yl - L y^, = L z^, (r=l, 2, . . . , n) , the eigenvalue problem (32) 

reduces to the standard form 


(o 2 y' = K'y', u) 2 z' ** K'z' 
r~r £r r~r -r 


r = 1, 2, 


where K'=L _1 KL T is a real symmetric positive definite matrix, in which 

L- T =(L-V. 

Denoting by v an arbitrary 2n-vector, Rayleigh's quotient associated 
with the eigenvalue problem (33) can be written in the form (Ref. 10) 

R(v) = v T K'v/v T v (34) 


Because K' is real and symmetric, it is well known that Rayleigh's 
quotient has a stationary value in the neighborhood of an eigenvalue, so 
that a stationarity principle exists also for gyroscopic systems. 

Now, we are in the position to discuss the truncation problem. The 
2n x 2n matrix K' was obtained as a result of representing the structure 
by an n-degree-of-freedora system. This representation is tantamount to 
the imposition of a given number of constraints on the original structure. 


For example, the constraints imposed on the system by the first of Eqs. 

(26) are T) c n +i = N +2 = **' “ 0. Truncating the series for u c by 
C C 
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assuming that i)_ „ - 0, we obtain a matrix K" obtained from K' by 

«» B C 

deleting two rows and the corresponding two columns. If the eigenvalues 

2 2 2 2 2 
gj of K' are such that a), < o) < ... < id and the eigenvalues 8 of K" 
r 2 2 1 — * j — n r 

are such that 0^ £ 0 2 — ’ • • £ ® n -l’ tben by fbe * nc l usj - on pri '-^tple 

(Ref. 3), we have 

2 2 2 2 
W 1 1 6 i 1 w 2 - 6 2 - ' 


. 2 „ fl 2 . 2 
< to .<8 , < w 

— n-1 — n-1 — n 


(35) 


Note that the fact that the eigenvalues of K* and K" have multiplicity 

two Is automatically taken into account in inequalities (35) . On the 

other hand, by relaxing one constraint, i.e. by adding one term to the 

series for u , we obtain a (2n+2)*(2n+2) matrix K"' which is obtained by 

2 

adding two rows and columns to K 1 . The eigenvalues of K ,M are such 
that 


2 2 2 2 
o A <. <a 2 <_ (d 2 < 


2 2 2 
<a <w <a: 1 
— n — n — n+1 


(36) 


The above developments permit us to conclude that the system 
natural frequencies tend to decrease monotonically with each additional 
degree of freedom . At the same time there is a new frequency added 
which is higher than any of the previous ones. Inequalities (35) and 
(36) have been derived by substracting and adding one term to the series 
for u^, i.e. the first of Eqs. (26). It is obvious that similar in- 
equalities will be obtained by subtracting and adding one term to the 
series for u^ or to the series for u^. Hence, it is possible to assess 
the effects on the system natural frequencies of truncating each sub- 
structure individually. 

It should be pointed out that the existence of an inclusion prin- 
ciple for gyroscopic systems is predicated on the existence of a sta- 
tionarity principle for such systems. 

7. Selection of Admissible Functions 

The question remains as to how to select the admissible functions. 
The first thing that comes to mind is to select them as the eigenfunc- 
tions of the various substructures. This approach has two drawbacks: 

1) the solution of the substructure eigenvalue problem may not be so 
easy to obtain and 2) an eigenvalue problem for the substructure may 
not be readily definable because of various rotational coordinates 
coupling the problem to that for a different substructure. Whereas 
substructure eigenfunctions (if they can be defined) would make a suit- 
able set of admissible functions, fortunately, there are many other sets 
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of functions that can yield equally good results, provided the functions 
satisfy certain criteria, as discussed below. 

The representation of substructure elastic motion in terms of 
admissible functions is a process closely related to the Rayleigh-Ritz 
method. Hence, it is natural to expect that the criteria for the 
selection of admissible functions; for the Rayleigh-Ritz problem associated 
with a particular single structure should be equally valid here, where 
the structure is now a substructure of a larger assembled structure. 
Sufficient criteria for each closed domain D are (see Ref. ll): 

1. Any finite number of admissible functions must be linearly 
independent. 

2. The set of admissible functions must be complete in the 
energy space of the domain D. 

3. The set of admissible functions should be strongly minimal 
in the energy space H A of the domain 5. 

The second of these criteria requires that we be able to identify 
a set of admissible functions that is complete in the energy space H . 

To this end, the following statement is useful: If A and B are positive 

definite operators and the spaces and Hg contain the same elements, 
then any set that is complete in H is complete in H . Therefore, by 
letting B have a simple form for which a complete set is easily found, 
we can use that set of admissible functions for the operator A. This 
justifies the earlier statement that it is not necessary to use substruc- 
ture eigenfunctions. Indeed, for an Involved operator A it should be 
possible to find a very simple operator B such that and Hg contain 
the same elements and a complete set of admissible functions for is 
found readily without any of the complexities inherent in working with 
the operator A. 

The third of the above criteria requires explanation. To this end. 


we consider an infinite system of admissible functions tp^j • 
in some space H. Then, the system is said to be strongly 


minimal in H if the smallest eigenvalue of the matrix 


R n * [( *i' Vh 1 


i.i 


n 

1 


(37) 


is bounded below by a positive number that is independent of n. An 


orthonormal system is strongly minimal by definition because the 


eigenvalues of R^ are equal to one for all n. Moreover, if operators 


16 



A and B are self-adjoint and positive definite In some Hilbert space H 
and every member of H. belongs to H_, then any set of admissible func- 
tlons in which is 
H. . Hence 

A 



8. Discrete Substructures 

In the preceding discussion, there was an Implicit assumption that 
the substructures were continuous. The same approach, however. Is valid 
for discrete substructures, as the dynamic characteristics are similar 
and do not depend on the mathematical model. Clearly, for discrete 
substructures no discretization is necessary but more often than not 
truncation is. For example, some discrete substructures can be repre- 
sented by mathematical models possessing hundreds or even thousands of 
degrees of freedom. Yet it is not feasible to work with systems with 
such a large number of degrees of freedom and truncation is an absolute 
necessity. This can be done by representing the motion of the substruc- 
ture by an adequate number of "modes". Typically, one would derive the 
eigenvalue problem for the substructure and attempt a solution, at least 
for a limited number of lower modes, and use these to represent the 
substructure. Of course, the question remains as to the "boundary 
conditions" to be imposed on the substructure. In the case of the 
central substructure one can regard it as entirely unrestrained. On the 
other hand, the choice of coordinates dictates that the substructure A 
be regarded as cantilevered at the origin A of axes x^y^z^. 

If the substructure possesses a very large number of degrees of 
freedom, then the solution of the eigenvalue problem can become not only 
time consuming, but the results inaccurate. Then, the question arises 
whether it is possible to avoid a solution of the eigenvalue problem 
altogether. In this regard, the authors of this report are trying to 
develop the concept of "admissible vectors", which are large-dimensional 
vectors analogous to the admissible functions discussed earlier. Ques- 
tions of definitions and convergence still remain, but the concept is 
sufficiently promising to warrant <1 sustained research effort. 

9 . Ulus trs tive Examp le 

Let us consider a specific structure consisting of a central sub- 
structure with a four-bladed rotor attached to it via a flexible 


, a system that is orthonormal in H_ is strongly minimal in 


strongly minimal in Hg Is also strongly minimal in 
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shaft. We shall take the central substructure In the fora of a beam and 
assume that the rotor and shaft rotate uniformly relative to the beam 
about an axis perpendicular to It (see Fig. 2). The central beam, the 
rotor blades, and the shaft are assumed to be one-dimensional members. 
The shaft can be considered as substructure A and the entire rotor as 
substructure B so that the shaft and rotor can be considered as » type 
AB chain of substructures. Moreover, assuming that the rotor has more 
than two blades and that the rotor blades are identical, so that the 
rotor is symmetric when undeformed, the resulting equations of motion 
for the entire structure can be reduced to a set of constant-coefficient 
equations representing a gyroscopic system. 

The actual derivation of the system kinetic energy is extremely 
involved and will not be presented here. Instead, we shall merely 
discuss the qutnMties describing the specific substructures which are 
needed in the kinematlcal considerations. A method for obtaining an 
explicit expression for the kinetic energy of the entire structure is 
discussed in the following sections. 

First, to describe the motion of the central beam, let us consider 


the set of axes x c yQZ c with the origin C at the mass center of the un- 
deformed beam and with the x axis along the length of the beam. We 

v 

shall take axes y and z r to be perpendicular to the beam so as to form 

u L* 

a triad in which the z axis is parallel to the shaft when the structure 
is undeformed. The orientation of x y z r relative to the inertial axes 

u V w 


XYZ is described by three rotations: 


6 C1 about x c> about y c> and 


0^2 about z c » in that order. 


so that the matrix T„_ of direction cosines 

CO 



where c0 cl and s0 cl denote cos 0 cl and sin © cl , respectively, etc. 
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Moreover, the radius vector from 0 to C, w^,, in terns of components 
along axes XYZ is 

*OC * ^ W 0C W 0C W 0C ^ ^ 39i 

x y z 


and the angular velocity of axes x c y c z c in terms of components along 
these axes can be shown to be 


c6 C2 C0 C3 
?C “ " c0 C2 s6 C3 


'®C3* 


The position vector of any mass point in the one-dimensional beam is 
simply 

r c - (x c 0 0> T (41) 

where x^ is the position of the point relative to C along the x^ axis 
and the elastic displacement of that point measured relative to x^y^z^, 

L li v 


u c = u c <x c , t) “ (0 v c (x c ,t) w c (x c ,t)} (42) 

Using Eqs. (38)-(42), we can obtain an explicit form of Eq. (1) for the 
absolute velocity of a mass point in the central beam. Such an explicit 
expression is not presented here for brevity. 

Next, let us consider the shaft, which is a type A substructure. 

We shall take the shaft to be connected to the central beam at point A 
a distance r from C along the x_ axis when the central beam is unde- 
formed, positive in the -x^, direction, and consider the axes x ca^ca z CA 
with origin at A and axis along the length of the shaft when it is 
undeformed. Assuming that axes ^^CA^CA are P ara H e ^- to a xes x^y^z^, 
when the central beam is undeformed, the matrix of direction cosines L QA 
between these axes due to geometry alone is the identity matrix. Hence, 
substituting Eq. (42) into Eq. (4) and evaluating the result at x c ■ 

- rCA, we can write the matrix of direction cosines L„. L„. between axes 
X CA^CA Z CA an< * a f ter deformation of the central beam as 


C 7 C C 


l ea l ga 


V c'^* r CA' t> "c'< C CA ,1: * 


- V c’<- r CA- t) 

-• U c' ( - r CA- C) 


where primes denote differentiations with respect to Similarly, we 

Can write the angular velocity vector 6 of axes x-.y_.z_. relative 

~LA LA LA LA 
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to Vc z c 
!ca' 


due to deformation in the form 

{0 -w c '(- r CA ,t) 


( 44 ) 


Assuming that in equilibrium the shaft axes x^y^z^ rotate relative to 
X CA y CA Z CA w * t *' t * ie un i^ orin angular velocity a about the z^ A axis, and 
that there is no oscillation of X A Y A Z A relative to ^^CA^A* t * ie 
angular velocity of x^y^z^ relative to x c A Yc A z c A * s 

“A “ <1, A L A L EA L GA T C0®3 (45) 

where e» * tO 0 1} and the matrix of direction cosines L. between 
^ J A 

Wa and x CA y CA z CA 18 8l ”" >ly 


“cos UKt 
A 

-sin to fc 
A 

- 0 


sin ui.t 
A 

t cos UK t 
A 


0 

0 (46) 

0 1 . 

In addition, the position vector of any mass point in the one-dimensional 
shaft can be written as 

0 z.r (47) 


8 a - t0 


Z A>‘ 


and the elastic displacement of that point measured relative to axes 
x^y^z^ has the form 

?A ■ ?A ( V 0 * tu A ( V E) v A (z A' t) ° ,T <48) 

Using Eqs. (38)-(48), we can obtain an explicit form of Eq. (7) for the 
absolute velocity of a mass point in the shaft. Again, such an expres- 
sion is omitted here for brevity. 

Finally, let us consider the entire rotor, whose geometric center 
is' connected to the shaft at point B, a distance r A g = L g from A along 

the z axis when the shaft is undeformed. Considering a set of axes 
A 

X AB Jf AB Z AB origin at B and parallel to axes X A Y A Z A when the shaft is 

undeformed, the matrix of direction cosines between axes x.y.z. and 

GB AAA 

x.„y.„z. D due to geometry is the identity matrix. Substituting Eq. (48) 
into Eq. (9), differentiating the result with respect to time, and 
evaluating the result at z *= L yields the angular velocity vector 

A S 

8^ of axes x A gy A g z A jj relative to axes x A y A z j due to the bending defor- 
mation. If, in addition, we allow the one- dimensional shaft to have 
torsion- 1 deformation, then the angular velocity vector 8^ has the 
explicit form 

?AB - { - i A ,(l .- t> “ A ' <I V t) *<V t))T <«> 

where primes here denote differentiation with respect to z^ and 
<f>(L g ,t) is the shaft torsional displacement at point B. Assuming 
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that each rotor blade is cantilevered to the shaft at point B, i.e., 
that axes x B y B z B do not rotate relative to axes x ^ B yAB z AB’ t * ,e an g u ^ ar 
velocity of x B y B z B relative to x ^ B y^g z ^ B is z eto and the matrix of 
direction cosines Lg between these two sets of axes is the identity 
matrix. Hence, the matrix of direction cosines Tg A between axes x B y B z B 
and x A y A z A can be shown to be 


BA EB 


^(L s ,t) 


■ i.. - -• »a..t) 


-“A ,<I v t) 


u A ’a s .t) v A '(L s ,t) 


To describe the position of a mass point in the rotor, we shall 
take one undeformed rotor blade to lie along the Xg axis, another to lie 
along the y fi axis, the third to lie along the -Xg axis, and the fourth 
to lie along the -y n axis. Hence, we consider four one-dimensional 

D 

members which are symmetric about the poinr B and have the domains of 
extension 0 <_ x g <_ Lg, y g = 0; 0 <_ y B _< Lg, Xg « 0; -Lg <_ x g <_ 0, y g = 

0; -Lg <_ yg 0, Xg = 0. We can write the position of a mass point in 
the domains of extension 0 < x < L , y ** 0 and -L < x <0, y = 0 

D — D 5> D D O 

as 

r g = {x B 0 0} T (51) 

and the position of a mass point in the domains of extension 0 _< Xg 

l b* y B = 0 and “ l b - x r - y B = 0 as 

r B = {0 y B 0} T (52) 

The elastic deformations of a mass point in each domain of extension 
measured relative to axes x B y B z B can be written as 

“B1 = "Bl (x B* t) = { ° V B1 (X B’ C> W Bl ( V t))T ’ 0 -1 X B - V y B = 0 


“ B 2 = ?B2 (y B ?t) = {-v B 2 (y B* C) ° W B2 (y B« t)} ’ 0 - y B - L B’ X B = ° 


B - B’ B 


'ui ( x n> t) (0 -v B3^ x g ,t ^ — X R — yR ® 


,B3 _B3 ' B 


B3 B' 


B - B 


?B4 = "B4 (y B» 0 = {v B4 (y B lt:) ° W B4 (y B’ t)} * _L B - y B^°> X B = ° 


where v n ,(5,t) 0 < x, < L n ; i = 1, 2, 3, 4 are the elastic displacements 

of each member in the x B y B z B plane and Wg^(?,t) 0 <_ X, <_ Lg; 1-1, 2, 3, 

4 are the elastic displacements of each member perpendicular to the 

x y z plane. Note that z, represents a local spatial coordinate for 
ii b D 
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each domain of extension. Using Eqs. (38)- (53), we can write an ex* 
pllcit form of Eq. (10) for the absolute velocity of any mass point in 
the rotor. As before, the resulting lengthy expression is omitted for 
brevity. 

Turning our attention to the system potential energy, we note that 
for the particular example under consideration axes x^y^Zg and 
do not oscillate relative to axes and X CA^CA Z CA* res P ect d ve ly » 

so that the 1 unction U in Eq. (14) does not appear in the potential 
energy. Moreover, assuming that the uniform angular velocity w is 

Pi 

large compared to all other angular velocity components, we can write 

Eq. (14) in the explicit form 

.2 


.. i f rT /* rt-'W 




dz. 


J 2,2 2 , 

+ “a <V l b ' 5 > 




d 5 


(54) 


where EI„, El,, and EI„ are bending stiffnesses of the central beam, the 
shaft, and each rotor blade, respectively, GJ^ is the torsional stiffness 
of the shaft, and p b is the mass per unit length of each rotor blade. 

It remains to specify the geometric boundary conditions for each substruc- 
ture. To this and, we note that the mass center of the undeforraed 
central beam always coincides with the origin of axes X ^V C Z C SD that 
u^(0,t) =0. In addition, we shall consider the shaft to be cantilevered 

to the central beam so that u,(0,t) = 1 = 0 and 4>(0,t) = 0. 

_A oza i ZA~" 

Finally, since each rotor blade is cantilevered to the shaft, we have 

u Bi ^u,t; l ?= o ~ 

Equations (42), (48), and (53) give the elastic displacement of a 
mass point in each continuous substructure in terms of continuous 
functions of the spatial position and time. As indicated in Sec. 5, 
it is necessary to eliminate the spatial dependence. To this end. 
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let us Introduce the specific discretization scheme for the example 
structure 



' 0 

0 

0 
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0 

0 

*B1* C * 

^B2^- 

(W 0 ) 

where 

[0 1 - denotes the 

3x4 null 

matrix 

. Note 

that we have 

taken n 

C = 4 ’ 


n, = 4, and n = 16 in Eqs. (54). In addition, to discretize the tor- 

A J3 

sional displacement (J>(z A ,t), we write 

* ( V t) - G 1 (z A )n <i)1 (t) + 0 2 (z A )n^ 2 (t) 

= {0 l (z A } 0 2 (z A )}{n n (t) n 4.2 (t) >T = (56) 

so that the entire structure is represented by 6+4+4+2+16 =32 
degrees of freedom. Of course the space dependent functions appearing 
in Eqs. (55) must satisfy the geometric boundary conditions mentioned 
previously. 
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10. G eneral Concepts for Symbolic Manipulation 

The major obstacle encountered in the derivation of explicit equa- 
tions of motion is the explicit evaluation of the system kinetic energy. 
Indeed, introducing Eqs. (1), (7), and (10) into Eq. (UK we obtain 

T = 2 m c?oc Yoc + 2 ?c J c?c + 2 J » “c“: dm c + ”oc T co"c 


- "oS’wcBc - h\ + I "A“cI»CA + 7 ?I J A?A 

1/ •T* • T T • *T T — — T 

+ ■7 / u u dm + w T u. - - H fi 

2 */ m A ..A..A A _CA AC.A „CA AC A_A _A^A 


+ f Was + \ + \fv^ "l" 3 dn B + i 


T T* 

T u 
AB BA_B 


T T— -T 

- _*AB T BA~ r B?B - » 


m c dra c * i c = J 
m c m c 


f (t c + Q c )dm c 


'~C=f "r dm C ’ J c = f (f c + Q c )T(E c + a c )dm c 

m m 


-f (t c + V’V'Y 

w m 


and similar definitions hold for ra. , ijl , r A , r n , u , u , J , J , H. , 
and Hg. Upon contemplating the form of Eq. (57), one concludes that 
evaluating the system kinetic energy explicitly involves the calculation 
of lengthy matrix products. Moreover, these matrix products involve 
quantities of different degrees. Generally, higher degrees imply less 
significance and, in fact, many higher-degree terms can be ignored. 
Specifically, for the linearized problem discussed herein, terms that 
are higher than second degree in the generalized coordinates and general- 
ized velocities are ignored. The complexity associated with the algebraic 
multiplication- of many matrices and the high probability of human error 
in performing many such multiplications by hand makes it highly desir- 
able to computerize these algebraic operations, a procedure whereby the 
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required algebraic operations can be Implemented on a computer Is 
outlined in this section. 

Let us first point out that obtaining Eqs. (30) in explicit form is 

X 

tantamount to determining the coefficient matrices M, F -F, and K^-K^ 
which in turn are determined by specific knowledge of the matrices M, F, 
and K^. The matrices M, F. and K^, can be obtained for a given 
structure by evaluating Eq. (57) explicitly in terms of the discretized 
coordinates, expanding the result in a Taylor series about the equil- 
ibrium configuration, applying a coordinate transformation to eliminate 
explicit time dependence, and then simply identifying the matrix ele- 
ments as the coefficients of the various products of generalized co- 
ordinates and generalized velocities. Of course, the matrix Ky is 
specified for a given structure because it is assembled from the sub- 
structure stiffnesses. 

The ideas used in implementing the required algebraic manipulation 
on a computer are best introduced via an explicit example. Let us con- 
sider the product 

0,5 m A"0C T C0 ^ T C0”0C “ (r CA + U CA^C + ?CA^ (59) 

1 * T* 

which appears in the expansion of y m ^ w cA W CA’ w ^ ere tlie l atter was en ~ 
countered in Eq. (57). The explicit calculation of this product is done 

in two steps. First, the resulting matrix expression is obtained by 

» • 

considering the vectors and u^ and the matrices T CQ , f and 

u CA symbolically, and then, knowing the elements of each matrix or 
vector explicitly for a specific structure, the matrix expression 
previously obtained is evaluated in terms of the explicit matrix ele- 
ments. By separating the calculation into two steps, it is possible to 
ignore certain matrix products based on the minimum degree associated 
with each matrix in the product without having to evaluate the matrix 
product itself explicitly. In addition, sometimes the result of 
evaluating a matrix product is known apriori. In this case, each 
occurrence of the matrix product in a matrix expression can be recog- 
nized immediately and the suitable substitution made. Hence, simpli- 
fication of the matrix expression is possible before it is evaluated 
explicitly, ultimately reducing the computation time required to obtain 
the explicit result. 
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For the moment, let us consider only those developments necessary 
for obtaining the matrix result. To calculate this product on the 
computer, we shall associate numbers with the different groups of 
characters that represent quantities to be manipulated algebraically, 
i.e., m A , w^, « c , u CA , T CQ , r CA , and u^. These groups of characters 
constitute what will be referred to as symbols . Note that some symbols 
such as m A are scalar symbols and others such as w oc » T cq» etc, » aTe 
matrix symbols . The association of numbers with symbols allows us to 
substitute manipulation of numbers or manipulation of symbols and is 
accomplished by the formation of a scalar symbol table and a matrix 
symbol table , Tables I and IX. The scalar symbol table contains two 
entries per line. These entries give the character representation of a 
matrix symbol, the weight assigned to the matrix symbol, the dimensions 
N and M of the N*M matrix symbol, and a pointer which locates the alge- 
braic representations of the individual matrix elements. This fifth 
column of the matrix symbol table containing the pointers will be dis- 
cussed later. The weight of a scalar (matrix) symbol is determined by 
the analyst according to his knowledge of the scalar (matrix) symbol's 
magnitude or his desire, to retain the symbol in the equations of motion. 
A high numerical value of the weight implies a less significant symbol. 

Note that here we assigned to the scalar m a weight of zero, to the 3*1 

A 

m • 

matrices w^, and u CA a weight of unity, and to the 3x3 matrices 

Tco» r^, aTU * we i8^ ts zero, zero, and unity, respectively. The 

number associated with a particular scalar or matrix symbol is the line 
number in the scalar or matrix symbol table. In this case, one is 
associated with m A , 51 with w Q £, 52 with etc. Note that the first 
line in the matrix symbol table is numbered consecutively after the last 
line in the scalar symbol table. 

Examining Eq. (59), we see that it is necessary to form algebraic 

•IT 

products of symbols, for example 0.5 m ^ w oc^'co* en< ^» we define 

terms . A term consists of a signed numerical coefficient , a pattern 
consisting of the numbers associated with each symbol in the product, 
and a weight representing the sum of the weights of each individual 
symbol appearing in the term. Note that the order of occurrence in a 
pattern of scalar symbols is arbitrary while the order of occurrence of 
matrix symbols must be preserved. Moreover, it is necessary to dis- 
tinguish the tranpose of a matrix symbol from the matrix symbol itself. 
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This is accomplished by inserting in the pattern at the appropriate 

place the negative of the number associated with the particular matrix 

symbol. In view of the number-symbol associations and the weights 

• T T 

of Tables I and II, we represent 0.5 m A w oc T CQ as a terra having 3 coef- 
ficient of 0.5, a pattern of 1, -51, -54 and a weight of 1. All terms, 
i.e. all coefficients, weights and patterns are stored in numbered 
storage stacks. The coefficient and weight of any term are always 
single numbers and are stored side-by-side in the coefficients and 
weights free storage stack , exhibited in the form of Table III. On the 
other hand, the pattern of a terra may differ from the pattern of other 
terms and must be able to represent the product of any number of sym- 
bols. Because of the different lengths of different patterns, all 

patterns are stored in the separate patterns free storage stack , labeled 

• T T 

as Table IV. Note that the coefficient and weight of 0.5 m A w 0 c T C0 are 
stored in line 3 of Table III and the pattern is stored in lines 5, 6 
and 7 of Table IV. 

It is also necessary to form algebraic sums of terms, such as 

• • 

T C0 W 0C _ ^CA^C~^CA n c +U CA* w h ic h we shall call series . We consider first 
only matrix series . Each matrik series is given a distinct series name . 

As examples, we shall call the matrix series consisting of the single 
term 0.5 by the name X, the matrix series T co *o C -^a?C~ S CaV?CA 

by the name Y, and the matrix series resulting from the product of X and 
Y by the name Z. A matrix series is described by a sequence of terms 
with coefficients and weights stored sequentially in the coefficients 
and weights free storage stack and with patterns stored sequentially in 
the patterns free storage stack. The summation of terms in the sequence 
is understood. Hence, the matrix series Y is described by the coefficients 
and weights stored in lines 5, 6, 7 and 8 of Table III and the patterns 
stored in lines 11—12, 14-15, 17-18 and 20 of Table XV. To distinguish 
sequences of terms forming matrix series, each matrix series is assigned 
a line in the matrix series definition table . Table V. Each line of 
this table contains seven entries giving: (1) the matrix series name, 

(2) the line number of the coefficient and weight in the coefficients 
ana weights ?.ree storage stack of the first term in the matrix series, 

(3) the number of terms in the matrix series, (4) the line number of the 
beginning of the pattern in the patterns free storage stack of the first 
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term in the matrix series, (5) and (6) the dimensions N and M of the N*M 

t 

matrix series, and (7) a pointer which locates the algebraic representa- 
tions of the individual matrix elements which would result if the matrix 
series were to be expanded explicitly. The seventh column in the 
matrix series definition table containing the pointers will be discussed 
later. Assigning line 3 in the matrix series definition table to the 
matrix series Y, the third line of Table V contains the entries 5 and 11 
giving the storage locations of the matrix series, the entry 4 indicating 
that there are four terms in the series, and the entries 3 and 1 indicat- 
ing that evaluation of the matrix series results in a three by one 
matrix. 

Let us now consider the multiplication of two matrix series, 
namely, the multiplication of Y by X, which can be performed term by 
term. The product of two terms yields a new term. If the total weight 

of the new term, given by adding up the weights of the two terms in the 

product, is greater than a specified value, for example 2, then the new 
term is deleted. Otherwise, the coefficient of the new term is the 
product of coefficients and the pattern of the new term is the conca- 
tenation of the patterns of the two terms in the product, where the 
concatenation preserves the order of occurrence of matrix symbols. 
Multiplying the four terms in Y by uhe single term in X, the first two 
resulting new terms each has a total weight of 2, the third new term has 
a total weight of 3 and is deleted, and the fourth new term has a total 
weight of 2. The first new term has a coefficient of 0.5 and a pattern 
of 1, -51, -54, 54, 51. The second new term has a coefficient of -0.5 
and a pattern of 1, -51, -54, 55, 52. The third new term is deleted and 
the fourth new term has a coefficient of 0.5 and a pattern of 1, -51, - 

54, 53. Assigning line 6 of Table V to the matrix series Z, the result- 

ing three terms retained in the product of Y by X are stored according 
to the information in line 6 of Table V. 

Turning our attention to the additional developments necessary for 
evaluating the matrix series obtained above explicitly (in terms of the 
algebraic representations of the individual matrix elements) and for 
storing the result, let us consider only the last term in the matrix 
series Z, namely 
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( 60 ) 


°- 5 ”a”oc t co“ca 

and discuss a matrix series Z1 formed from this single term and stored 
according to the information in line 7 of Table V. To evaluate Zl, Eq. 
(60), explicitly, the 3 x 1 matrices w^ and u^ and the 3*3 matrix T^ 
must be known explicitly. To this end, we shall use 


"oCy\ ’ ,CA 


rCAl\l +<,) CA2\2 

*CAl 1 * 1 C3 +<f ’CA2 n C4 


’1 0 O' 
0 10 
0 0 1 


where the elements of these three matrix symbols represent scalar 
series and < 0 ^, ^ n cl » n c2 , ^ 3 , i QCx , w QCy and w QCz are 

scalar symbols with weights as defined in lines 2-10 of Table I. 

Scalar series are algebraic sums of terms with the pattern of each term 
containing only scalar symbols. Like matrix series, a scalar series is 
also described by a sequence of terms with coefficients and weights 
stored sequentially in the coefficients and weights free storage stack 
and with patterns stored sequentially in the patterns free storage 
stack. Again, the summation of terms in the sequence is understood. To 
distinguish sequences of terms forming scalar series, each scalar series 
is assigned a line in the scalar series definition table , Table VI. 

Each line of this table contains three entries giving the line number of 
the coefficient and weight in the coefficients and weights free storage 
stack of the first term in the scalar series, the number of terms in the 
series, and the line number of the beginning of the pattern in the 
patterns free storage stack of the first term in the scalar series. 

Note that no series name is associated with a scalar series. This is 


because scalar series are always algebraic representations of particular 
matrix elements and, as such, they are directly associated with either a 
matrix symbol or a matrix series which already has a series name. 
Considering the 3x1 matrix symbol w QC and its explicit definition given 
in Eq. (61), the scalar series representing the (1, 1), (2, 1) and (3, 

1 ) matrix elements are stored according to the information found in 
lines 1, 2 and 3 of the scalar series definition table. Moreover, the 

association between these three lines in Table VI and the matrix symbol 

• 

Wq£ is accomplished by entering 1, giving the location in Table VI of 
the (1, 1 ) matrix element, in column 5 of line 51 of the matrix symbol 
table. Table II. 
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If the matrix symbol in question has more than one column, then the 
information defining the scalar series representations of the matrix 
elements is stored in column-order in Table VI, i.e., one column after 
another. For example, the scalar series representing the (1, 1), (2, 

1), (3, 1), (1, 2), (2, 2), (3, 2), (1, 3), (2, 3) and (3, 3) matrix 
elements of the 3 X 3 matrix symbol are stored according to the infor- 
mation found in lines 4, 5, 6, 7, 8, 9, 10, 11 and 12, respectively, of 
Table VI. Note the 4 appearing in column 5, line 54, of Table II which 
locates the (1, 1) matrix element of T . 

*» • 

Now that the elements of the matrices w^, u^ and T^ are scalar 
series defined explicitly and associated with the appropriate matrix 
symbols through column five of the matrix symbol table, the explicit 
evaluation of the matrix series Z1 can be performed according to the 
usual rules governing matrix multiplication. The evaluation requires 
simply the multiplication and addition of scalar series. The multipli- 
cation of scalar series is exactly analogous to the multiplication of 
matrix series discussed earlier and is performed term by term, where 
again any new terms with a total weight greater than 2 are deleted. The 
result of the explicit evaluation of Z1 is a scalar series representing 
the (1, 1) element of the resulting 1x1 matrix. This scalar series is 
stored according to the information found in line 16 of the scalar 
series definition table and is associated with the matrix series Z1 by 
entering 16 in column 7, line 7, of Table V. 

11. The Computer Program 

The ideas outlined in the preceeding section are appealing because 
of their simplicity and can be programmed easily in FORTRAN. Indeed, a 
FORTRAN computer program implementing these ideas has been written. 

The program is organized into subroutines, with each subroutine per- 
forming a different operation required in producing the matrices M, F 
and K^, explicitly (such as multiplication of matrix series, multiplica- 
tion and addition of scalar series, etc.). However, in order to produce 
these matrices explicitly for a given structure, it is necessary to for- 
mulate an algorithm which executes the appropriate subroutine at the 
appropriate time. It is expected that the algorithm will be slightly 
different for each specific structure. In order to facilitate the 
analyst's job in formulating and changing the algorithm for a structure, 
a simple interpreter has also been written (in FORTRAN). The interpreter 
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accepts as input any of a set of twelve simple Instructions and, 
depending on the instruction, causes a particular subroutine or compli- 
cated sequence of subroutines to be executed. Two of the twelve in- 
structions are used to declare a group of characters to be either a 
scalar symbol with a weight or a matrix symbol with a weight and dimen- 
sions. Four instructions are used to input and output either a matrix 
series or the scalar series which define a matrix symbol or matrix 
series explicitly. One instruction is used to set the maximum weight 
for retaining terms. Another instruction is used to substitute a 
matrix series for every occurrence of a particular matrix symbol in 
another matrix series. Still another instruction is used to cause a 
matrix series to be evaluated in terms of the explicit matrix elements 
for each matrix symbol. Two instructions are used to define the vector 
of generalized coordinates and the vector of generalized velocities in 
terms of the individual scalar symbols representing each generalized 
coordinate and generalized velocity. And finally, one instruction is 
used to generate the coefficient matrices M, F and K^, from the scalar 
series representing the system kinetic energy in explicit form. 

12. Example Application 

Let us consider the specific structure discussed in Sec. 9 con- 
sisting of a central beam with a four-bladed rotor attached to it via a 
flexible shaft. The specific quantities needed in the kinematical 
procedure for this structure are presented in Sec. 9. In this section, 
however, for brevity of presentation, we shall represent each continuous 
coordinate by only one degree of freedom. This is tantamount to sub- 
stituting ^ C 2^ x c^ :=< ^A2^ Z A^ =< * > B2^^ =0 in Eqs * and 0 2 ^ z a^ =O in Eq * 

(56)* Hence, the entire structure is represented by 6+2+2+1+8=19 degrees 

of freedom and the 19-dimensional configuration vector is given by 
q = * W OCx’ W OCy’ W OCz’ 9 C1’ 0 C2’ e c3’ n Cl’ n C3‘ n Al’ n A3’ n <j>l’ 

T 

n Bl’ n B3’ n B5’ n B7’ n B9* n Bll’ n B13’ n B15 } (62 

To derive explicit equations for small motions about equilibrium, the 
developments of Sec. S are used to calculate the mass integrals, Eqs. 
(58). These integrals are presented for each substructure in Appendix 
A. Using the computer program and these mass integrals in conjunction 
with the kinematical relationships given in Sec. 8, the system kixietic 
energy, Eq. (57), has been evaluated explicitly and found to depend on 
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time. In order to eliminate the explicit time dependence, we substitute 
the coordinate transformations 


n Al " n Al COS V + n A3 SlnU A t 


n A3 = ' n Al sln V + n A3 C ° 8W A t 


(63a) 


n Bl = 

n* — 

r B1 

'$5 + 

n* 9 coscu A t -f- 

n Bl3 S±nW A t 

n B3 

n* - 

n B3 

n B7 + 

n Bll C ° S V 4 

’ n B15 S±nW A t 

n B5 = 

n* + 
n Bl 

^85 

n B9 SinW A t + 

n^cos^t 

n B7 = 

n* + 
] B3 

n* - 

n B7 

n Bll S±n V 4 

’ n B15 COS V 

n B9 “ 

n* - 
r B1 

n B5 ~ 

n*9 c °sm A t - 

n|i 3 sinm A t 

n Bll = 

n"k — 

n B3 

n B7 ~ 

r 'Bll COS< °A t 

" n B 15 S±n V 

n B13 = 

n Bl + 

' n B5 + 

n* 9 sinm A t - 

n B13 C ° S V 

n E15 

n B3 + 

^7 + 

n Bll Sin V 

- n* 15 cosm A t 


(63b) 


into the mass integrals of Appendix A. The resulting mass integrals are 
presented in Appendix B. Note that Eqs. (63b) represent a transfor- 
mation to multiblade coordinates (see Ref. 5). The computer program is 
now used to evaluate the system kinetic energy explicitly in terms of 
the transformed coordinates and then to identify the coefficient matrices 
M, F and K^,. The constant matrices M, F and K^, along with the constant 
matrix Ky obtained from Eq. (54) by considering Eqs. (55) and (56) in 
conjunction with Eqs. (63) are presented element by element in Appendix 
C. 

To obtain Eq. (30) from the matrices M, F, K^, and found In 
Appendix C, it remains to consider the elimination of ignorable coor- 
dinates as discussed in Sec. 4. For the particular structure under con- 
sideration, the vector w^ and the rotational coordinate are ignorable 
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coordinates. Regarding the conservation of the generalized momenta 
associated with these coordinates as constraints on the system, ve can 
use these constraints to eliminate and 9^ from the formulation and, 
as a result, obtain the equatons of motion for the system in the form of 
Eq. (30) . The constraint equations in question are presented in Appen- 
dix D. 


13. Concluding Remarks 

The purpose of this investigation is to develop a general approach 
to the dynamic synthesis of substructures consisting of an assemblage of 
point-connected flexible substructures. At least one of the substruc- 
tures undergoes uniform rotation in the equilibrium state, so that the 
system is gyroscopic. The mathematical model is intended to permit the 
determination c*f the modal characteristics of helicopters. 

The approach is based on the substructure synthesis concept, where- 

* 

by the elastic motion of a flexible substructure is represented by a 
suitable set of "modes". These modes need not represent natural modes 
of the substructures (corresponding to given boundary conditions) but 
can consist of a set of so-called admissible functions. A stationarity 
principle and an ensuing inclusion principle for gyroscopic systems 
permit the assessment of the effect of truncating the set of admissible 
functions on the system natural frequencies. Some rational criteria for 
the selection of the admissible functions are provided. 

The method is applied to a simplified model of a helicopter, con- 
sisting of a beam simulating the airframe, a transmission shaft, and 
four elastic blades simulating the rotor. Use of the method of multi- 
blade coordinates permits the reduction of the system of equations with 
periodic coefficients to one with constant coefficients, where the 
latter is of gyroscopic type. 

The enormous amount of algebraic operations involved in the deriva- 
tion of the coefficient matrices makes automation a virtual necessity. 

To this end, a program for symbolic manipulation on a digital computer 
has been developed. The program has been used to derive the coefficient 
matrices for the specific example considered. 

This investigation represents a first step toward a rational ap- 
proach to a problem of extreme complexity. Many aspects of the problem 


still require answers. One of these problems is the concept of admis- 
sible vectors. This would permit the representation of discrete substruc- 
tures by a set of vectors which do not necessarily represent natural 
modes of the substructure. Another problem requiring attention is that 
of multipoint connections for flexible gyroscopic systems. Although 
this problem may not appear to be so critical in a helicopter, the 
capability of simulating fixed wings should be beneficial. 
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FIGURE I. THE MATHEMATICAL MODEL 
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32766 denotes the end of a term 

32767 denoted the end of a series 

Table VI - Pattern Free Storage Stack 
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Appendix A 

1. The mass integrals for substructure C 
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ii. The mass integrals for substructure A 
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M A = *A11^A1 + ^A3* 
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iii. The mass integrals for substructure B 
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Appendix B 

i. The mass integrals for substructure A 
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B A + *All (n A^ 2 “A t+n A3 s2 “A t+2, 'Al ,l A3 S “A t “A t) 
» 2 Al (n Al 5 "A t -"A3 C “A t) 


•^ Al < ri Al C “ A t +,| A 3 S V ) 1 

♦ 2 Al ( ’ , L 9 “ A t - n A 3 ‘ : " A t > 

- , * 2 , * 2 . 

*All* n Al +n A3* 


(Bl) 


.*2 • *2 2 *2 *2 .**.** 

M A = 4 ’All^ n Al +n A3 +UJ A (n Al +n A3 ) + 2a) A^ n Al n A3 _n A3 T1 Al^ 


<B2) 


“A 


' * A1 [ n A1 C ‘ J A t+n A3 S ‘ 1) A t+w A ( - I1 Al SU A t+,1 A3 Ca) A t > ] ' 
^A1 t- ^Al Sa) A t+ ^ A t_l ° A (n A 1 C “ At+n A ' lSa,it J - 1 ' 


'A3 A A A1 A 
0 


A3 A 


(B3) 


r A 


1 T 
2 m A l S 


~ J m A L S 
0 


L*Al (n Al S “A t - n A3 C “A C) 


♦Al (n Al C “A e+n A3“A t) 


♦ ai ^ V ^ V ^ as ^ a 0 ' 

”^ Al ^ T 1 Al CtU A t +, 1 A 3 SU A t ^ 


0 


?A* 


'* 2 Al [ -"Al SU A t+ "A3 C “A t -“A ( ' 1 Il C “A ,:+n 23 S “A ,:)1 

5 zAl l ‘ B M“A t -"A3“A t ' to A ( "Al S “A ,: ‘ n A3 < ”A t)1 

W V ) ] 


<B4) 


(B5) 


where and were defined in Appendix A. 
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ii. The mass integrals for substructure B 
The elements of are: 

t , on ,,T I* Z 4.->r>* 2 + * 2 x9 * 2 x * 2 2 . *2 *2 2 -j. * 2 

J B11 * 2B B +2 *Bll [n Bl +2n B3 +ri B5 +2n B7 +n B9 C “A t+n BU +n B13 S V +n B15 
~ 2n Bl n B5' f2n BS ri B13 SU> A tClll A t3 
J B12 " ^zBl n B5 

^ ^ ^ A A 

J B13 ' - 2 *zM <ri Bn C “A t+n B15 S “4 t)+2 *Bll [ - n B3' , B9 S V +rl B3"B13 C “A t 

-"B7 n M a “A t+n B7C3 C V' n Bl"BU 9 "A t; ‘ n B5'’Bll S V +,| Bl ri B15 C “A t 
+ ' , B5 , ’b15 C “a' ;1 
J B21 “ 4 *zBl n B5 

*2 *2 *2 *2 *2 2 *2 *2 2 *2 * * 

J B22 ’ ZB B +2 *Bll l ’ 1 Bl +2n B3 +,1 B5 +2,l B7 + "B9 s “A t+n Bll +n B13 c “A t+,, B15 +2n Bl , ’B5 

- 2n B9 ,1 B13 s “A t “A 11 

J B23 ' 2 * 2 Bl <n BU S "A C -' , B15 C V )+2 *BU t -' 1 B3 n B9 C “A ,; ‘ n B3 , ’B13 S “A C+n B7 , 'B9 C “A t 

it it it it it it it it 

+n B7 n B13 SW A t-T1 Bl n Bll C(il A t:+n B5 n Bll CU A t '’ 0 Bl ri B15 SU) A t+n B5 n B15 Sai A t -* 


J"' 1 - = “ 2 ^zBl (n Bll C1J A t:+n B15 S(1, A t)+2 ^Bll [ " n B3 T1 B9 StU A t:+n B3 n B13 C(1J A t_n B7 ri B9 Sl0 A t 


B31 


it it it it it it it it it it 


* * 


* * 


* * 


J B32 = 24 ’zBl (n Bll St0 A t n Bl5 CU A t)+2<tl Bll^“ n B3 n B9 CllJ A t " n B3 n B13 SI11 A t+ri B7 n B9 C11) A t 


it it it it it it it it it it 

+T1 B7 n B13 S “A t: " n Bl n Bll Ca) A t+n B5 n Bll Gai A t ' i ' n Bl T1 B15 SW A t+n B5 n Bl5 S ' 1, A t;l 


*2 *9 *2 *2 , 

J B33 = 4B B +2lJ> Bll (2n Bl +2n B5 +n B9 +n B13 ) 


(B6) 
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Mg - 2$ 


Bll 


, »*2 «*2 **2 »*2 **2 **2 **2 »*2 2 . *2 *2 *2 *2 
^ 2n Bl +2n B3 +2n B5 +2n B7‘ fn B9 +n Bll +n Bl3 +n B15 +liJ A^ T1 B9 +n Bll +n B13 +n B15^ 


B5 
» * 


+2w A^B9 n Bl3”' 1 B13 r, B9 +n Bll n Bl5" n B15 ri Bll )1 


(B7) 


«B 


r 2 ^itn* 9 s V^B13 G V + ^ (ri B9 C V +n B13 SW A t) J * 
2 *Bl^B9 Ctl> A ,:+ ^BX3 aa, A t+tlJ A ( " n B9 SW A t:+n B13 CW A t: ^ 
4 ^B1^B3 


(B8) 


r B = 


4 ^Bl n B3 


^ 2i Bl (T, B9 ea, A t+n B13 S “A t) 


~ 4 *Bl n B3 


2 *B1< WV + WV 0 " 

2 ^Bl ( - n B9 Sa> A t+n B13 C V ) 


2 ‘ f> Bl^ n B9 Sa, A t_n B13 CU) A t ^ 


(B9) 


The components of the vector Hg are*. 

®B1 = 2 *zBl^ SaJ A t (n Bll +aJ A n B15 ) '* CU A t ( n Bl5“‘ U A n Bll ) J +2 *Bli* C * il A t ^ n Bl n Bll 


. .* * 


»A *•**•** •* * .* *.**.** 

■ n Bll n Br n B3 n B9 +n B9 n B3~ n B5 n Bll +n BU n B5 +n B7 n B9“ n B9 n B7 
+a] A (ri B3 n B13“ T1 Bl ri B15" n B7 n B13 +T1 B5 fl Bl5 ) ]+SW A t tn Bl n B15 _ ^B15 n Bl 

• A A «A A » A A • A A »A A »A A A A 

-n B3 n B13 +ri B13 n B3 -r, B5 ri B15 +n Bl5 ri B5 +n B7 T1 B13" ri B13 r, B7 +U) A (n Bl T1 Bll 

— A A A A A A 

n B3 n B9 +n B7 n B9“ n B5 n Bll^ 

®B2 " 2 + Z Bl IS “A t<; B15'VBll )+C “A t<; 'Bn' h “A n B15 ) 1+2 +Bll {c “A t[; 'Bl ,, B15 

.A A • A A .A A.AA • A A • A A .A A 

“ n B15 n Br n B3 n B13 +n B13 n B3 +n B5 n B15' n B15 n B5~ n B7 n B13 +n B13 n B7 
+U A < ‘ n Bl n Bll“ n B3 n B9 -n B7 n B9 +n B5 T, Bll^ ] +Sa * A t ^" n Bl n Bll +n Bli n Bl 

• A A .A A .A A .A A .A A .A A A A 

' n B9 n B3 +n B3 T1 B9 -ri B5 n Bll +n Bli n B5~ ri B7 r, B9 +n B9 ri B7 +W A ( ~ n B3 n B13 


+n Bl n B15 n B7 n B13 +T1 B5 ri B15 ) ^ 
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^B3 

where 


- 4 $ 


• A 

zBl n Bl 



*zBl’ 


and B_ were defined In Appendix A. 

OX b 


(BIO) 
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A ppendix C 

Throughout this appendix any element of a matrix which is not 
specified is equal to zero, 
i. The matrix M 

The elements M (r,s = 1,2,..., 19) of M have the following values 
rs 

M 11 ' M 22 ■ K 33 " b C + “a + 4 "b 


M 


15 51 


M r , = -M„, =* -M, „ = ~ m.L 0 + 


*24 


‘42 2 AS 


M 


44 ~A 
= M 


B. + 4nigLg + 2B 


B 


* m t x„ A + 


35 “53 “26 “62 “*A"CA ' 4m B x CA 

m 55 - B c + "a4 + b a + 4 “b <4 + 4 + zb b 


M 


M 


66 " B C + VCA + 4m B X CA + 4B B 

[ 46 " M 64 ‘ <5 “A + 4 "B ,1 ‘CA L S 

M 77 ” 4 C11 * <”a + 4 V*cia + 44 C1A B B 

H 67 * M 76 ' 4 xCl ‘'“A + 4 “b )x CA 4 C1A + 4 *C1A B B 

»«7 * M 74 * -<T »A + 4 V*CU 1 S 

M 27 = M 72 = M 38 = M 83 = ^Cl + (m A + 4m B )<J> ClA 

2 2 


m qo = 


,2 


88 3 ^Cll + (m A + 4m B H ClA + (B A + 4m B L S + 2 V*C1A 
M 58 = M 85 = -*XC1 + (m A + 4m B )X CA^ClA “ (B A + 4m B L S + 2 V*C1A 
M 18 = M 81 = “ * 2 m A + 4m B^ClA L S 
M 99 = ^All + 4m B^AlB + 2B B*A1B 
M 89 = M 98 = ~^zAl + 4 ”b L S*A 1B + 2B B^A1B^C1A 
M 59 = M 95 = " M 4,10 = ~ M 10,4 = *zAl + ^VaIB + 2B B <{ 'a1B 
M 91 3 M 19 = M 2,10 = M 10,2 = ^A1 + 4i Va 1B 
M 10,10 = ^All + 4m B^AlB + 2!(> A1B B B 
M 7,10 = M 10 , 7 3 ( ^A1 + 4m B^AlB^CIA 
M 6,10 ~ M 10,6 = ~ (m A^Al * 4 , Va 1B )x CA 

M ll,ll 3 40 aib b b 

M 7,ll = M ll,7 3 4 ^C1A G A1B B B 


48 



ORIGINAL PAGE IS 
OF POOR QUALIW 

**6,11 “ **11,6 " 4G A1B B B 


**12,12 * M 13,13 " **14,14 * **15,15 “ 4 *B11 

**7,12 “ M 12,7 “ 4 *€1A*zB1 

** 6, 12 * ** 12,6 = 4 ^ zB 1 

**8,13 = **13,8 ’ 4 *C1A*B1 

**5,13 ^ **13,5 = 4x CA^B1 

**3,13 * **13,3 = 4 *B1 

**16,16 = **17,17 “ **18,18 * **19,19 " 2 *B11 
**10 , 16 = **16,10 = "**9,18 = “**18,9 “ 2<,i A1B*B1 


**7,16 = **16,7 = 2 *C1A*B1 
**8,17 = **17,8 " 2<J> C1A*zB1 


**9,17 = **17,9 = **10,19 “ **19,10 * " 2 *A1B*zB1 
**8,18 = **18,8 = 2 *C1A L S*B1 
**6,16 = **16,6 = “ 2x CA^B1 
4,16 16,4 5,18 18,5 S B1 


M = M = - M 
2,16 n 16,2 1,18 


“ **18,1 ~ 2 *B1 


**5,17 **17,5 " **4,19 “ **19,4 _ “ 2 *zB1 


where we have used the definitions 


*C1A = 5 *Cl (x C ) l x = -x 

C C.A 


^AIB 3 *Al (z A^z. « L_ 
A S 


. d<l ci (x c ) | 

*C1A ~ 3x_ I 

C X C ' * X CA 

. . d V (z A ) | 

A1B ^ ■* - 4 


0 A1B 0 A1 (z A^z * L 
A S 

All other quantities are defined in Appendix A or the text. 
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ii, The matrix F/oia 


The elements F (r,s » 1,2,..., 19) of F/w have the fr 

ITS A 


values : 


F 41 " F 52 - ( 2 m A + 4m B )L S 


F 45 “ B A + 4m B L S + 2B B 
.2 


F 48 “ - (B A + 4m B L S + 2 V*C1A 

F 49 = F 5,10 = ^zAl + 4m B L S*AlB + 2 *A1B B B 

F 4,17 = F 5,19 = “ 2 *zB1 

F 4,18 = “ F 5,16 = “ 2 ^B1 L S 

F 54 ' ' B A - 4 "B L S + 2B B 

F 56 ■ - ( T “A + 4 ”b ) ’‘ca l s 

F 57 ’ ( l “A + 4m B ) +ClA L S 

F 82 = " ^2 m A + 4m B^ 4 ClA L S 

F 84 ■ <B A + 4 "B L S - 2 V*C1A 

F 86 = m A + 4m B )x CA^ClH L S 

F 87 = ^ 2 m A + 4m B^ClA^ClA L S 

**.10 = “ ^zAl + 4m B L S*AlB + ^AIbV^CIA 

F 8,16 = " 29 C1A L S^B1 

F S, 19 = 2<? C1A^zB1 


92 


F 10,l “ 4m B^AlB 


F 94 = F 10,5 “ 4m B^AlB L 8 + 44, A1B B B 
F 96 = 4m B X CA^AlB 
F 97 ~ ~ 4m B^ClA^AlB 
F 9,10 = ” F 10,9 = 


4m B*AlB 


F 9,16 = F 10,18 = “ 2<j> AlB*Bl 
F 10,8 = “ 4< ' m B L S <5) AlB AIbV^CIA 
F 17,4 = F 19 ,5 “ ~ 4 ^zBl 

F 16,18 = F 17, 19 = - F 18 , 16 = " F 19,17 = “ 2 ^B11 
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F 17 t 10 F 19,9 * «W«1 

F 19,8 “ **ClA*zBl 


2 

iii. The matrix Kt/^a 

2 

The elements (r,s = 1,2,..., 19) of have the following 

values : 

444 = 4ss = b a + 4m B L S " 2l3 B 

K T4 , 10 = *4:10,4 = " *4:59 “ " 4:95 = 4m B^AlB L b + 4 b b*aib 

44,19 “ 4.19,4 = " 45,17 “ " 4l7,5 = " 4 *zBl 

458 = 485 = " <B A + 4m B L S “ 2B B**C1A 

488 = < B A + 4m B L S + 2 V*C1A 

4s 9 = 498 = 4m B^ClA L S l{, AlB + 4 *C1A*A1B B B 

48.17 = 4l7 ,8 = " 4,}> C1 a4b1 

499 = 4l0,l0 = 4m B^AlB 

4] 2,12 = 4l4, 14 ~4l6,16 = 4l8,18 = “*B11 

417 . 17 = 4l9, 19 = 2 4ll 


iv. The matrix Ky 

The elements (r,s = 1,2,..., 19) of have the following 

values; 


477 ” 488 




dx. 


, . , _ r\t (AAAX 

«V99 - “VIO.IO J EI A^ 2 ) 

A 

S/ll,!! - f G \ dz A ) Z A 


dz. 


4l2,12 4l3 , 13 " 4l4,14 ~ 4l5,15 _ 4l6,16 2 4l7,17 “ 2 4l8,18 

Jo t £I B 


2 4l9, 


/d ? 2 2, / d *Bl <C) \ 

) + "a4 (l b - < > \~d ?— ) 


19 


d^ B1 (c)\ 2 


Id? 
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Appendix D 

The constraint equations are 

<"c + m A + % ) “0Cx + "A + 4m B )L S ( “c2 + VW ' <1 "a + '"V*C1aV'c3 
+ <*A1 + SW"Al * 2 *Bl\l3 + 4 VA1BVAB ” X 1 * D1) 

(« c + » A + 4 -B)”oCy - <7 “A + 4m B )L S ( ®Cl ” “a 6 C2> ' <"a + 4m B )x CA°C3 

+ + VclA + ‘VclA^C! - "A + 4 WclAVc3 


+ ( *A1 + 4n B*AlB );i A3 - ‘VaIbVaI + 2 *BlV * X 2 <D2 > 


(m c + m A + AW* + <m A + 4B B )x CA e C2 + ( +Cl + m A+ClA + 4m B*ClA ) n C3 
+ 4 ^Bl n B3 ~ X 3 


- (» A + 4® B > x CA”oCy + ( 2 “A + 4 "b )x CA L S W C1 * “a S C 2 ) + (B C + Vc* 

+ 4 VCA + 4 V°C3 + (5 XC1 - n A x CA*ClA - 4m B X CA 4 ClA + ^'cuV'cl 
* ( I "A + 4 V’ t CA*ClA L s“A n C3 " <m A*Al + 4 Va1B >X CA S A3 
+ 4, VcA*AlB“A n Al + 4 °A1B B bV + 4? zBl"Bl ‘ x CA*BlV * X 4 <(D4) 


where all quantities are defined in Appendix A, Appendix C and the text. 
The constants (i = 1,2, 3, 4) are determined by means of the initial 
conditions. Quite often they are taken as zero. 
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